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We present a detailed comparison of computational efficiency and precision for several free energy 
difference (Af) methods. The analysis includes both equilibrium and non-equilibrium approaches, 
and distinguishes between uni-directional and bi-directional methodologies. We are primarily inter- 
ested in comparing two recently proposed approaches, adaptive integration and single-ensemble path 
sampling, to more established methodologies. As test cases, we study relative solvation free ener- 
gies, of large changes to the size or charge of a Lennard- Jones particle in explicit water. The results 
show that, for the systems used in this study, both adaptive integration and path sampling offer 
unique advantages over the more traditional approaches. Specifically, adaptive integration is found 
to provide very precise long-simulation AF estimates as compared to other methods used in this 
report, while also offering rapid estimation of A_F. The results demonstrate that the adaptive inte- 
gration approach is the best overall method for the systems studied here. The single-ensemble path 
sampling approach is found to be superior to ordinary Jarzynski averaging for the uni-directional, 
"fast-growth" non-equilibrium case. Closer examination of the path sampling approach on a two- 
dimensional system suggests it may be the overall method of choice when conformational sampling 
barriers are high. However, it appears that the free energy landscapes for the systems used in this 
study have rather modest configurational sampling barriers. 
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I. INTRODUCTION 

Free energy difference (AF) calculations are useful for 
a wide variety of applications, including drug design-'^'S., 
solubility of small molecules^''*, and protein/ligand bind- 
ing afBnities^iSii. Due to the high computational cost of 
AF calculations, it is of interest to carefully compare the 
efficiencies of the various approaches. 

We are particularly interested in assessing recently 
proposed methods*'^ in comparison to established tech- 
niques. Thus, the purpose of this study is to provide 
a careful comparison of the efficiency and precision of 
several A_F methods. We seek to answer two important 
questions: (i) Given a fixed amount of computational 
time (10^ dynamics steps, in this study), which method 
estimates the correct value of AF with the greatest preci- 
sion? (ii) Which AF approach can obtain a "reasonable" 
estimate of AF in the least amount of computational 
time? 

Free energy difference methods can be classified as 
either equilibrium or non-equilibrium. Equilibrium ap- 
proaches include multi-stage free energy perturbation^", 
thermodynamic integrationiii^, Bennett analysis^^'^** 
and weighted histogram analysia^^. The common theme 
in these approaches is that sufficiently long equilibrium 
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simulations are performed at each intermediate stage of 
the free energy calculation. Equilibrium methods are 
in wide use and are known to provide accurate results; 
however, the computational cost can be large due the 
simulation time needed to attain equilibrium at each in- 
termediate stage. A host of non-equilibrium methods 
have recently been applied to various molecular systems, 
largely due to Jarzynski's remarkable equalitjiiSiii. Non- 
equilibrium methods have the potential to provide very 
rapid estimates of AF, but can suffer from significant 
bias^&^. 

In this report we present results using both equilib- 
rium and non-equilibrium approaches — as well as uni- 
directional and bi-directional methodology. Specifically, 
we compare: (i) adaptive integration^; (ii) thermo- 
dynamic integration^^; (iii) single-ensemble path sam- 
pling of non-equilibrium work values using Jarzynski's 
uni-directional averaging^; (iv) single-ensemble path 
sampling using Bennett's bi-directional formalism; (v) 
Jarzynski averaging of non-equilibrium work valuesiSiSi; 

(vi) Bennett analysis of non-equilibrium work valuesi^i^S; 

(vii) equilibrium Bennett analysisi^ii^; and (viii) multi- 
stage free energy perturbatior*iS. We also compare the 
free energy profiles, which determines the potential of 
mean force, for adaptive integration and thermodynamic 
integration. 

Generally, one is interested in the free energy difference 
(AF — Fi—Fq) between two states or systems of interest 
denoted by potential energy functions Uo{x) and Ui{x), 
where x is the full set of configurational coordinates. AF 
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can be written in terms of the partition functions for each 
state 



AF = -ksTln 



' z[Ui{x)\ 



(1) 



where fc^ is the Boltzmann constant, T is the system tem- 
perature, and Z[J7(x)] = / c?xexp[— [/(a;)/fcBT]. Because 
the overlap between the configurations in C/g a-nd Ui may 
be poor, a "path" connecting C/q and C/i is typically cre- 
ated. In our notation, the path will be parameterized 
using the variable A, with < A < 1. 



II. EQUILIBRIUM FREE ENERGY 
CALCULATION 

Equilibrium free energy methodologies share the com- 
mon strategy of generating equilibrium ensembles of con- 
figurations at multiple values of the scaling parameter 
A. In the current study we investigate thermodynamic 
integrationli, adaptive integration^, multi-stage free en- 
ergy perturbationifi, and multi-stage equilibrium Bennett 
analysisi^. We performed separate equilibrium simula- 
tions at successive values of A, and then estimated AF 
using free energy perturbation, Bennett averaging, and 
thermodynamic integration on the resulting ensemble of 
configurations (detailed in Sec. lIV|l . 



A. Thermodynamic integration 

Thermodynamic integration (TI) is probably the most 
common fully equilibrium AF approach. In TI, equilib- 
rium simulations are performed at multiple values of A. 
Then, AF is found by approximating the integralii, 



AF = 



A=0 



dUxjx) 
dX 



(2) 



where the functional form for U\{x) depends upon the 
scaling methodology and will be discussed in detail in 
Sec. II VI The notation (...}a indicates an ensemble aver- 
age at a particular value of A. In addition to the possibil- 
ity of inadequate equilibrium sampling at each A value, 
error arises in TI from the fact that only a finite number 
of A values can be simulated, and thus the integral must 
be approximated by a sumM. Thermodynamic integra- 
tion can provide very accurate AF calculations, but can 
also be computationally expensive due to the equilibrium 
sampling required at each A valu8Si2iiSS42&. 



B. Adaptive integration 

The adaptive integration method (AIM), detailed in 
Ref. 0, seeks to estimate the same integral as that of 
TI; namely Eq. ((2Jl (see also discussions in Refs. I?^l2^ 



However, in addition to fixed- A equilibrium 
sampling, the AIM approach uses a Metropolis Monte 
Carlo procedure to generate equilibrium ensembles for 
the set of A values. The A-sampling is done by attempting 
Monte Carlo moves that change the value of A during the 
simulation. The probability of accepting a change from 
the old value Ao to a new value A„ is 



-^acc(Ao ^ Aj^ 



1.0, e 



, (3) 



where /3 = l/ksT and 5F{Xi) is the current running free 
energy estimate obtained by numerically approximating 
the integral 



SF{X,) = 



A=0 



(4) 



Between attempted Monte Carlo moves in A, any canoni- 
cal sampling scheme (e.g., molecular dynamics, Langevin 
dynamics, Monte Carlo) can be used to propagate the 
system at fixed A. In this report, Langevin dynamics is 
used to sample configurations, and Monte Carlo moves 
in A are attempted after every time step. 

It is important to note that, due to the use of the 
running estimate 5F in Eq. (PJ , the AIM method satisfies 
detailed balance only asymptotically. In other words, 
once the AF estimate fully converges, the value of SF is 
correct, and detailed balance is satisfiedSi^. 

AIM is related to parallel tempering simulationSi, and 
has the associated advantage: equilibrium sampling of 
conformational space at one A value can assist sampling 
at other A values due to the frequent A moves. This 
is reminiscent of "A dynamics" simulationS^iSS, but con- 
trasts with TI where only a single starting configuration 
is passed between A values. 

An additional advantage of AIM over the other meth- 
ods detailed in this report is that there is a simple, built- 
in, reliable, convergence criterion. Specifically, one can 
keep track of the population (number of simulation snap- 
shots) at each value of A. When the estimate for AF has 
converged, the population will be approximately uniform 
across all values of A. If the population is not approxi- 
mately uniform, then the simulation must be continued. 



C. Free energy perturbation 

In the free energy perturbation approach, one performs 
independent equilibrium simulations at each A value (like 
TI), then uses exponential averaging to determine the 
free energy difference between neighboring A values^^ — 
these differences are then summed to obtain the total free 
energy difference. AF can be approximated for a path 
containing n A- values (including A = and A = 1) using 
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the "forward" estimate (FEPF) 

ri-l 

or the "reverse" estimate (FEPR) 

n-l 



i=0 



A primary limitation of free energy perturbation is that 
the spacing between A values must be small enough that 
there is sufficient overlap between all pairs (Ai, A^+i) of 
configuration spaces. 



D. Equilibrium Bennett estimation 

It is also possible to use Bennett's method to combine 
the information normally used for forward and reverse 
free energy perturbation. In this approach, one computes 
the free energy difference between successive A values SFi 
according to 



1 + f,P{u^.+iiSi)-Ux,{S.}-SF,) 



(7) 



Ai+1 



Then the sum of these 6Fi is the total free energy 
differencei^; 



AF = ^ 6F,. 



(8) 



i=0 



Studies have shown that using the Bennett method to 
evaluate free energy data is the most efficient manner to 
utilize two equilibrium ensemblesi^i^. 



III. NON-EQUILIBRIUM FREE ENERGY 
ESTIMATION 



In non-equilibrium free energy approaches, the system 
is forced to switch to subsequent A values, whether or 
not equilibrium has been reached at the current A value. 
In this way, non-equilibrium paths are generated that 
connect Uq and Ui. In the current study we use uni- 
directional Jarzynski averaging^ and bi-directional Ben- 
nett averaging of Jarzynski-style work valuesl^., as well 
as uni-directional^ and bi-directional averaging of path 
sampled work values. 



A. Jarzynski averaging 

For the Jarzynski method^^, one considers non- 
equilibrium paths that alternate between increments in A 
and "traditional" dynamics (e.g., Monte Carlo or molec- 
ular dynamics) in x at fixed A values. Thus, a path with 
n A-steps is given by 

Z„ = |(Ao = 0,xq), (Ai,xo), (Ai,xi), (A2,xi), 

(A2,X2), (A„_i,f„_i), (A„ = l,x„„i)|, (9) 

where it should be noted that increments (steps) from 
Xi to Ai+i are performed at a fixed conformation Xi, and 
the initial xq is drawn from the canonical Uq distribution. 
For simplicity, Eq. (|5J| shows only a single dynamics step 
performed at each fixed Xi, from Xi-i to Xi] However, 
multiple steps may be implemented, as below fSec. IVTl. 
A "forward" work value is thus given by 



n— 1 



(10) 



By generating multiple paths (and thus work values) 
it is possible to estimate AF via Jarzynski's equality^^ 



AF = -kBT\ii{e-l^^*-)^. 



(11) 



where the (...)g represents an average over forward work 
values Wf generated by starting the system at Uq and 
ending at C/i . A similar expression can be written for the 
situation when work values are generated by switching 
from Ui to Uq. This approach is "uni-directional" since 
only work values from cither forward or reverse data are 
used. 

Perhaps the most remarkable aspect of Eq. Hll|) is that 
it is valid for arbitrary switching speed. However, in prac- 
tice, the AF estimates are very sensitive to the distribu- 
tion of work values, which in turn is largely dependent on 
the switching speed. If the distribution of work values is 
non-Gaussian and the width is large {ayy ^ ksT), then 
the AF estimate can be heavily biasediSiiSiSSiS. Consis- 
tent with results in this report (Sec.0), other efficiency 
studies^^'^° have suggested that the optimal efficiency for 
uni-directional Jarzynski averaging is when the switching 
speed is slow enough that aw ~ 1 fcsF. 



B. Bennett averaging of Jarzynski work values 

Due to the bias introduced in using uni-directional 
Jarzynski averaging, it is useful to consider a method 
where both forward and reverse work values are uti- 
lized. It has been shown that the most efficient use of 
bi-directional data is via Bennett's methodii^S, 



E 

Nt 



j{-n+Wt-\F) 



E 



f}[-r)+W,+AF) 



in) 
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where rj = ksT In {N{/N^) allows for differing number 
of forward (Nf) and reverse {N,-) work values. Equation 
(I12|l must be solved iteratively since AF appears in the 
sum on both sides of the equation. 



C. Single-ensemble path sampling 

Single-ensemble path sampling (SEPS) is a non- 
equilibrium approach that seeks to generate "important" 
paths more frequently^sMiSSi^SiSLS. The method uses im- 
portance sampling to generate paths (and thus work val- 
ues) according to an arbitrary distribution D, here chosen 



D(Z„) = Q(Z„) e 



(13) 



where Q(Z„) is proportional to the probability of occur- 
rence of an ordinary Jarzynski path, and is given below. 
With this choice of D the free energy is estimated via 
(compare to Refs. 34 35 36 37 38) 



r 



(14) 



where the is a reminder that the work values used in 
the sum must be generated according to the distribution 
in Eq. (|13() . Since forward work values, Wf are utilized 
in Eq. H14|l . the paths must start in Uo and end in J7i. A 
similar expression can be written for reverse work values 

To generate work values according to the distribution 
D, path sampling must be used9i34i35i36,37,38,39,40,4i^ 

path sampling, entire paths are generated and then ac- 
cepted or rejected according to a suitable Monte Carlo 
criteria. In general, the probability of accepting a trial 
path with n A-steps (Z^j with work value W') that was 
generated from an existing path (Z„ with work value W) 
is given by 



pZ„- 
^ acc 



QiK) Pi< 



,cn 



•Z„ ^-^f3W' 



Q(Z„) Pt 



gen 



(15) 



where P^^^ is the conditional probability of generating 
a trial path Y from existing path X . 

For this study, we generate trial paths by randomly 
choosing a "shoot" point Ag along an existing path (com- 
pare to Refs. ,40.42.43,) . Then, Langevin dynamics is 
used to propagate the system from As — > (backward 
segment), followed by As ^ 1 (forward segment). Be- 
fore running the backward segment, the velocities at the 
shoot point must be reversed and then ordinary Langevin 
dynamics are used to propagate the system^°. Once the 
trial path is complete, all the velocities for the backward 
segment are reversed. Since the stochastic Langevin al- 
gorithm is employed in the simulation, it is not necessary 
to perturb the configurational coordinates at the shoot 
point to obtain a trial path that differs from the existing 
path. 



The above recipe for generating trial paths leads to the 
following statistical weights for the existing Q(Z„) and 
trial (3(ZJJ paths 



g(Z„) = e-'5^o(Jo) -Q ^ 



1=0 
n-l 



(16) 



1=0 



where p{xi —>■ x^+i) is the the transition probability for 
taking a dynamics step from configuration Xi to x^+i— . 
We have assumed for simplicity that only one dynamics 
step is taken at each value of A; however, the approach 
allows for multiple steps. The corresponding generating 
probabilities for the existing and trial paths are given by 



z„^z' 



Pchoosc Ppcrturb P{^i ~* ^i+l) JJ^ Pl^^i+l ^ ^j)? 



i=0 



pz;,^z„ 



n-l 



Pchooso Ppcrturb H P^^' ^ H ^ ^0,(17) 

i—s i—0 

where p{xi+i —> Xi) is the transition probability of tak- 
ing a backward step from x^+i to Xi. The "bar" notation 
is a reminder that the velocities are reversed for these 
segments. The probability of choosing a particular shoot 
point As is denoted by Pchoose, and the probability of 
a particular perturbation to the configurational coordi- 
nates at the shoot point is given by Pporturb- 

Since we have chosen not to perturb the configura- 
tional coordinates at the shoot point, and any value of A 
along the path is equally likely to be chosen as the shoot 

point, then Pporturb = Pporturb and Pchoosc = p'choosc- 
addition, since the transition probabilities obey detailed 
balance and preserve the canonical distribution then^^ 



p{Xi+i X 



) ^p{x, ^ f,+i) e"^(^^»+i("'^-^^»+i'"'+^0(I8) 



Inserting Eqs. (HHJl, and lO into Eq. lO gives the 
acceptance criterion for trial paths (compare to Eq. (45) 
in Ref.lSF 



+z' 



l,e 



-I3(^SW-5W' + ^{W'-W)) 



(19) 



where SW is defined as the work accumulated up to the 
shoot point for the existing path 



s — 1 

'5^ = Il[t^A.+i(x.)-C/A.(x.) 



(20) 



and SW' is the equivalent quantity for the trial path. 
Note that Eq. H19|) is independent of the details of the 
fixed- A dynamics. 
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To clarify ambiguities in our original presentation of 
the SEPS approach^, we also give details for applying 
it using overdamped Langevin dynamics (i.e., Brownian 
dynamics). In Ref. '5, backward segments were gener- 
ated using ordinary dynamics with negative forces, i.e., 
to be very clear, the force was taken to be identical to 
the physical force, but opposite in sign. Thus, the tran- 
sition probabilities for forward and backward steps are 
approximately equal 



p{Xi+l Xi) K.p{xi Xi+i). 

(Brownian dynamics) 



(21) 



Equality occurs when the forces at Xi and Xi+i are iden- 
tical. The acceptance criterion becomes 



r'a.cc = mm 



1 ^-PihiW -W) + Uo{S'q)-Uo{3o)) 



(Brownian dynamics) (22) 

Therefore, the criticism raised in a recent paper— is in- 
correct. 



D. Bennett averaging of path sampled work values 

The use of bi-directional data is worth considering 
for the SEPS method, just as it was for ordinary non- 
equilibrium Jarzynski work values. Generalizing Ben- 
nett's method to include the work values sampled from 
D gives 



E 



D (. + hPWt 



N{ 1 + e 



N, 



Nr 1 + e' 



P(-V+W, + AF) 



Nr 



(23) 



Thus, to obtain a Bennett-averaged estimate for AF, the 
path sampling algorithm is applied to generate an ensem- 
ble of paths going from to Ui (Wf , forward) and also 
for Ui to Uq (Wr, reverse). Then, Eq. (^5)1 is applied to 
the data. 



of the system was maintained at 300.0 K using Langevin 
dynamics with a friction coefficient of 5.0 ps~^. RAT- 
TLE was used to constrain all hydrogens to their ideal 
lengths'*®, allowing a 2.0 fs time step. A cutoff of 12.465 
A was chosen for electrostatic and van-der-Waals inter- 
actions with a smoothing function implemented from 
10.465 to 12.465 A. It is expected that the use of cutoffs 
will introduce systematic errors into the Ai^ calculation, 
however, in this report we are only interested in compar- 
ing Ai^ methodologies — we do not compare our results 
to experimental data. 

For the first test neutral Lennard- Jones par- 

ticle was "grown" from 2.126452 A to 6.715999 A. The 
sizes were chosen to be that of lithium and cesium from 
the OPLS-AA forcefield^^. In the second test case, the 
Lennard-Jones particle remains at a fixed size of 2.126452 
A, but the charge is changed from -e/2 to -f e/2. For each 
test case, and each AF method, the system was initially 
equilibrated for 100 ps (5 x 10"* dynamics steps). The 
initial equilibration is not included in the total compu- 
tational time listed in the results, however, since every 
method was given identical initial equilibration times, the 
efficiency analysis is fair. 

The A-scaling (i.e., the form of the hybrid potential 
U\) used for all AF methods in this study was chosen 
to be the default implementation within the TINKER 
packageiS=. If a particle's charge is varied from go to gi, 
the hybrid potential is simply the regular potential en- 
ergy calculated using a hybrid charge of 



q\ = Ml + (1 - ^)9o. 



(24) 



Similarly, if a particle has a change in the van der Waals 
parameters r, e the hybrid parameters are given by 



r\ = Ari + (1 - A)ro, 
fx = Aei -I- (1 - A)eo. 



(25) 



The free energy slope as a function of A for both the 
growing and charging test cases arc shown in Figs. Q and 
121 The smoothness of both plots suggests that a more 
sophisticated A-scaling is not necessary for this study. 
If, for example, we had chosen to grow a particle from 
nothing, then it is likely that a different scaling would be 
needed (such as in Refs. 



IV. SIMULATION DETAILS 

To test the efficiency and precision of each method de- 
tailed above we use two relative solvation free energy 
calculations. One involves a large change in the van 
der Waals radius of a neutral particle in explicit solvent 
( "growing" ) , and the other is a large change in the charge 
of the particle while keeping the size fixed ( "charging" ) . 

The system used in both cases consists of a single 
Lennard-Jones particle in a 24.93 A box of 500 TIP3P wa- 
ter molecules. For all simulations, the molecular simula- 
tion package TINKER 4.2 was used^S.. The temperature 



A. Thermodynamic integration calculations 

For thermodynamic integration (TI), equilibrium 
simulations were performed at each value of A. 
An equal amount of simulation time was devoted 
to each of 21 equally spaced values of A = 
0.0, 0.05, 0.1, ...,0.9, 0.95, 1.0. Averages of the slope 
dF/dX — {dU/dX)\, shown in Figs. ^ *ind 13 were col- 
lected for each value of A. The first 50% of the slope 
data were discarded for equilibration. Finally, the data 
were used to estimate the integral in Eq. Q using the 
trapezoidal rule. Note that higher order integration 
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schemes were also attempted, but did not change the 
results, suggesting that the curves in Figs. ^ and 13 are 
smooth enough that high order integration schemes are 
not needed for this report. Also, the percentage of data 
that was discarded for equilibration was varied from 25- 
75% with no significant changes to the results. 



the last A = configuration — thus the A = equilibrium 
ensemble was generated "on the fly." 

Similarly, "reverse" non-equilibrium paths were gener- 
ated by starting each simulation from configurations in 
the Ui equilibrium ensemble and switching from A = 1 
to A = 0. 



B. Adaptive integration calculations 

Adaptive integration (AIM) results were obtained 
by collecting the slope of the free energy dF/dX = 
{dU/dX)x, by starting the simulation from an equili- 
brated configuration at A = and performing one dy- 
namics step. Immediately following the single step, a 
Monte Carlo move in A was attempted, which was ac- 
cepted with probability given by Eq. Q. The pattern 
of one dynamics step followed by one Monte Carlo trial 
move was repeated until a total of 10® dynamics steps 
(and thus 10® Monte Carlo attempts) had been per- 
formed. The same A values used in TI are also used for 
AIM, thus A = 0.0, 0.05, 0.1, 0.9, 0.95, 1.0 are the only 
allowed values. For this report Monte Carlo moves were 
attempted between neighboring values of A only, i.e., a 
move from A=0.35 to 0.4 or 0.3 may be attempted but 
not to 0.45. Also, all SF{Xi) values of Eq. were ini- 
tially set to zero. The estimate of the free energy was 
obtained by numerically approximating the integral in 
Eq. ^ using the trapezoidal rule. As with TI, higher 
order integration schemes did not change the results. 



C. Free energy perturbation and equilibrium 
Bennett calculations 

All free energy perturbation calculations (forward Eq. 
Q and reverse Eq. ©), and equilibrium Bennett com- 
putations (Eq. ||SJ)) were performed on the same set of 
configurations as for TI. Specifically, equilibrium simula- 
tions were performed at each of 21 equally spaced values 
of A = 0.0, 0.05, 0.1, 0.9, 0.95, 1.0, and the first 50% of 
the data were discarded for equilibration. 



D. Jarzynski estimate calculations 

Estimates of the free energy using the non-equilibrium 
work values were computed using Eq. (|ll|l for Jarzynski 
averaging, and Eq. (|12|l for Bennett averaging. "For- 
ward" non-equilibrium paths were generated by start- 
ing the simulation from an equilibrated configuration at 
A = 0, then incrementing the value of A, followed by an- 
other dynamics step, and so on until A = 1. Thus, only 
one dynamics step was performed at each value of A. The 
work value associated with the path was then computed 
using Eq. H10|) . Between each path, the system was sim- 
ulated for 100 dynamics steps at A = 0, starting with 



E. Single-ensemble path sampling calculations 

For the single-ensemble path sampling (SEPS) 
method, we first generated an initial path using stan- 
dard Jarzynski formalism. The only difference between 
the paths described above and the initial path for SEPS 
was that, due to the computer memory needed to store 
a path, the number of A-steps was limited to 500 for this 
study. In other words, if the desired path should con- 
tain around 2000 dynamics steps, the simulation would 
perform four dynamics steps at each A value giving a to- 
tal simulation time of 1996 dynamics steps for each path 
(note that simulation at A = 1 was not necessary). 

Once an initial path was generated as described above, 
a trial path was created by perturbing the old path as 
described in Sec. IIII CI Then, the new path was ac- 
cepted with probability given by Eq. H19II . Importantly, 
if the new path was rejected, then the old path was 
counted again in the path ensemble. Also, as with any 
Monte Carlo approach, an initial equilibration phase was 
needed. For this report, the necessary amount of equi- 
libration was determined by studying the dependence 
of the average free energy estimate, after 10® dynam- 
ics steps, from 16 independent trials, as a function of the 
number of paths that were discarded for equilibration. 
The optimal number of discarded paths was then chosen 
to be where the average free energy estimate no longer 
depends on the number of discarded paths. 



V. RESULTS AND DISCUSSION 

Using the simulation details described above, two rel- 
ative solvation free energy calculations were carried out 
in a box of 500 TIP3P water molecules. Each of the free 
energy methods described above were used to estimate 
AF. Specifically, we compare: 

• adaptive integration (AIM) using Eqs. ^ and ||2J); 

• thermodynamic integration (TI) using Eq. ((SJ; 

• uni-directional single-ensemble path sampling 
(SEPS) using Eq. 

• bi-directional single-ensemble path sampling with 
Bennett averaging (BSEPS) using Eq. 03 : 

• uni-directional Jarzynski averaging of work values 
(Jarz) using Eq. ifTTjl : 



7 



Steps AIM TI SEPS 


BSEPS Jarz BJarz Benn 


FEPF 


FEPR 


2E3 16.3(4.6) 16.5(6.1) — 


— — — 16.7(6.2) 


18.7(6.7) 


14.5(5.7) 


4E3 14.4(3.9) 13.2(4.4) — 


— — — 13.4(4.4) 


14.7(4.7) 


11.9(4.2) 


9E3 10.4(3.3) 11.2(3.6) — 


— 7.9(1.3) — 11.3(3.6) 


12.3(3.9) 


10.1(3.3) 


1.7E4 8.94(2.35) 9.7(2.46) — 


— 7.56(0.93) 7.53(1.13) 9.75(2.46) 10.48(2.70) 8.92(2.26) 


3.5E4 7.51(0.52) 8.32(1.35) — 


— 7.62(0.84) 7.47(0.71) 8.36(1.38) 


8.91(1.63) 


7.74(1.11) 


7E4 7.38(0.48) 7.89(1.17) — 


— 7.55(0.67) 7.38(0.59) 7.92(1.19) 


8.35(1.40) 


7.46(0.97) 


1.3E5 7.35(0.36) 7.18(0.65) 7.15(0.79) 


— 7.34(0.49) 7.36(0.38) 7.22(0.64) 


7.56(0.68) 


6.83(0.68) 


2.7E5 7.34(0.23) 7.19(0.22) 7.19(0.62) 6.95(0.56) 7.35(0.44) 7.28(0.24) 7.21(0.22) 


7.29(0.25) 


7.08(0.20) 


5.5E5 7.22(0.12) 7.18(0.11) 7.19(0.29) 7.12(0.46) 7.32(0.28) 7.23(0.20) 7.18(0.12) 


7.22(0.11) 


7.16(0.13) 


1E6 7.19(0.07) 7.26(0.18) 7.17(0.18) 7.23(0.20) 7.25(0.23) 7.22(0.14) 7.26(0.18) 


7.28(0.18) 


7.24(0.20) 



TABLE I: Free energy difference estimates obtained for changing the Lennard- Jones size of a neutral particle in a box of 
explicit water. Results are shown for various methods described in the text as a function of the number of dynamics steps used 
in the simulation. Table entries are the mean estimates from 16 independent simulations with the standard deviation shown in 
parentheses. For single-ensemble path sampling (SEPS and BSEPS) and Jarzynski methods (Jarz and BJarz), only the most 
efficient results are shown. The table shows that in the limit of long simulation times (10® dynamics steps) all methods produce 
average AF estimates that roughly agree. The table also shows that AIM provides the most precise long-simulation estimate. 



• bi-directional Bennett averaging of Jarzynski work 
values (BJarz) using Eq. (|12|l : 

• Equilibrium Bennett approach (Benn) using Eq. 
©; and 

• multi-stage free energy perturbation in the forward 
(FEPF) and reverse (FEPR) directions, using, re- 
spectively Eqs. 10 and ©. 



A. Growing a Lennard-Jones particle 

We first compute the free energy required to grow a 
neutral particle from 2.126452 A to 6.715999 A in 500 
TIP3P waters. 

Figure ^ shows the slope of the free energy {dF/ dX = 
{dU/dX)x) as a function of A for both TI and AIM after 
10^ Langevin dynamics steps. The figure suggests that 
AIM can more efficiently sample the profile. In AIM, 
configurations are not forced to remain at a particular A, 
but may switch to another value of A if it is favorable 
to do so. Such "cross-talk" is apparently the source of 
the smoother A-profile compared to TI. Table^shows AF 
estimates for the different approaches used in this report. 
Note that for all non-equilibrium approaches, only the 
most efhcient data are shown. For SEPS and BSEPS 
all paths were composed of 500 A-steps (restricted to 500 
due to computer memory) with 40 dynamics steps at each 
value of A. For Jarz and BJarz the paths were composed 
of 10 000 A-steps with one dynamics step at each value 
of A. For all of these non-equilibrium data, the standard 
deviation of the work values were aw ~ 0.8 kcal/mol ~ 
1.3 ksT, in agreement with previous studiesi^i^S. At 
least five different path lengths were attempted for each 
non-equilibrium method to determine the most efRcient. 

Table ^ demonstrates that, for long simulation times, 
all methods produce roughly the same average AF esti- 




0.2 0.4 0.6 Q.i 

Switching parameter, X 



FIG. 1: The slope of the free energy dF/dX as a function of 
A for changing the Lennard-Jones size of a neutral particle in 
a box of explicit water. Results for both TI and AIM meth- 
ods are shown for lO'' dynamics steps. The data show the 
averages (data points) and standard deviations (error bars) 
from 16 independent simulations for each method. The figure 
demonstrates that AIM has the ability to sample the A-path 
more efficiently, thus producing a much smoother and more 
precise profile compared to TI. Thus, AIM is preferred over 
TI for computing the potential of mean force for this sys- 
tem. In addition, the smoothness of the profile suggests that 
the switching function U\ of Eq. 1251 used in this report is 
adequate. 



mate. Also, the table clearly shows that, given 10^ dy- 
namics steps, AIM provides the most precise free energy 
estimates. 

Table m shows the approximate number of dynamics 
steps needed by each method to obtain a free energy esti- 
mate within a specific tolerance of Ai^iong sim (average of 
all estimates at 10^ dynamics steps). Note that the num- 
ber of dynamics steps needed for the SEPS and BSEPS 
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Method Within 1.0 kcal/mol Within 0.5 kcal/mol 



AIM 


23 000 


30 000 


TI 


89 000 


181 000 


SEPS 


140 000 


377 000 


BSEPS 


279 000 


444 000 


Jarz 


18 000 


127 000 


BJaiz 


26 000 


96 000 


Benn 


90 000 


180 000 


FEFF 


104 000 


191 000 


FEPR 


60 000 


184 000 



TABLE II: Number of dynamics steps necessary to be within 
a specified tolerance of the correct result Aflong sim ~ 7.23 
kcal/mol, average AF estimate at 10® dynamics steps for aU 
methods, for growing a Lennard- Jones particle in exphcit sol- 
vent. The first column is the method used to obtain the es- 
timate. The second column is the number of dynamics steps 
needed to estimate AF within 1.0 kcal/mol of AFiong sim with 
an uncertainty less than 1.0 kcal/mol. The third column is the 
number of dynamics steps needed to obtain an estimate within 
0.5 kcal/mol with an uncertainty less than 0.5 kcal/mol. 



methods are large due to the fact that whole paths must 
be discarded for equihbration of the path ensemble. For 
all methods except AIM, the table entries for Table ^ 
were estimated using linear interpolation of the data in 
Table From the data in Table ^ if the desired preci- 
sion is less than 1.0 kcal/mol, then AIM, Jarz and BJarz 
appear to be the best methods. However, if the desired 
precision is less than 0.5 kcal/mol, then AIM is the best 
choice. 

Tables |2 and m taken together, demonstrate the dif- 
ference between using equilibrium data in the "forward" 
(FEPF) and "reverse" (FEPR) directions. While, the 
results are similar for 10^ dynamics steps, it is clear 
that FEPR produces the desired results more rapidly 
than FEPF indicating that the configurational overlap 
is greater in the reverse direction. However, the FEPR 
data also tends to "overshoot" the correct value by a 
small margin which makes convergence of the FEPR es- 
timate difficult to judge. 

Thus, we conclude that, for growing a Lennard-Jones 
particle in explicit solvent, the preferred method depends 
upon the type of estimate one wishes to generate. If a 
very precise high-quality estimate is desired, then AIM 
is the best choice by a considerable margin. If a very 
rapid estimate of AF, with an uncertainty of less than 
1.0 kcal/mol, is desired, then then comparable results 
are seen using AIM, Jarz and BJarz methodologies. If 
the AF estimate is to be within 0.5 kcal/mol, then AIM 
is the best choice. 

Finally, if the desired result is the potential of mean 
force, then AIM will generate a much smoother curve 
than TI. 



1. Fast-growth uni- directional data 

We now consider non-equilibrium uni-directional fast- 
growth data, i.e., generated by switching the system 
rapidly from Uq (small particle) to J7i (large particle). 
Importantly, there will be an advantage to generating 
uni-directional data in some cases, since only the Uq equi- 
librium ensemble is needed to estimate AF. 

In contrast to the data shown in Tables HI and ITTl where 
the lengths of the non-equilibrium switching trajectories 
were pre-optimized, here we focus on the efficacy of the 
methods using non-optimal, rather fast switching. After 
all, when attempting a free energy computation on a new 
system, there is no way to know in advance the optimal 
path length (number of A-steps). Substantial optimiza- 
tion may be needed for both SEPS and Jarz methods to 
work efficiently. 

Here, we test the SEPS and Jarz methods using short 
paths with an equal number of dynamics steps. For 
SEPS, 500 A-steps with four dynamics steps at each value 
of A was used, producing a distribution of work values 
with aw = 2.1 kcal/mol. For Jarz, 2000 A-steps with 
one dynamics step at each value of A was used, produc- 
ing a distribution of work values with aw — 2.9 kcal/mol. 
Note that these paths are roughly ten times shorter than 
optimal and thus aw is 3-4 times larger than the optimal 
value of ~ ksT. 

Figure [21 gives a comparison between SEPS and Jarz 
methods for the fast-growth uni-directional paths just 
described. The upper figure (a) shows the average free 
energy estimates and standard deviations for both the 
SEPS and Jarz methods. The lower figure (b) gives the 
histogram of the work values for each method. Both fig- 
ures also show the "correct" value AFiongsim, generated 
from a very long simulation. The figures clearly demon- 
strate that, for fast-growth data, SEPS has the ability to 
"shift" the work values such that the AF value is near 
the center of the work value distribution — rather than 
in the tail of the distribution as with the Jarz method. 
Thus, the SEPS results converge more rapidly than Jarz 
to the correct value of AF. 

We suggest that the the SEPS method may find the 
greatest use for the ability to bias fast-growth work values 
to obtain the correct value of AF, as shown here. 

B. Charging a Lennard-Jones particle 

We next compute the free energy required to charge a 
particle from -e/2 to -|-e/2 in 500 TIP3P waters. 

Figure 13 shows the slope of the free energy (dF/dA = 
{dU/dX)\) as a function of A for both TI (green) and 
AIM (black) after 10^ Langevin dynamics steps. The 
data shown in the plot are the mean (data points) and 
standard deviation (errorbars) for 16 independent trials. 
While the errorbars are too small to resolve on the plot 
shown, the average uncertainty in the the slope for AIM 
is 0.38 kcal/mol and for TI is 1.05 kcal/mol, suggesting 
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FIG. 2; (a) "Fast-growth" uni-directional free energy difference estimates obtained for changing the Lennard- Jones size of a 
neutral particle in a box of explicit water. Results are shown for both SEPS and Jarz methods as a function of the number of 
dynamics steps used in the simulation. For both methods, fast-growth work values were generated by simulating roughly 2000 
dynamics steps per path, which is ten times shorter than optimal. The solid horizontal line represents the best estimate of the 
free energy difference AFiang aim based on averaging all results shown in Table |I] at 10^ dynamics steps. The averages (data 
points) and standard deviations (errorbars) are from 16 independent simulations, (b) Histograms of the work values used to 
generate the free energy estimates for both the SEPS and Jarz methods. The plots demonstrate the potential usefulness of 
using path sampling over regular Jarzynski averaging. Specifically, if the work values are fast-growth and uni-directional, then 
SEPS is able to bias the work values in such a way to improve the free energy estimate. Note that for all the SEPS data shown, 
the first 50 work values are thrown away for equilibration as, described in Sec. II V El 



Steps AIM TI SEPS 


BSEPS Jarz BJarz Benn FEPF FEPR 


2E3 8.5(5.5) 24.5(2.3) — 


— — — 24.4(2.3) 28.7(2.8) 20.0(2.1) 


4E3 9.7(6.6) 21.5(3.0) — 


— — — 21.4(3.1) 25.4(3.0) 17.7(3.1) 


9E3 14.6(11.4) 20.1(1.7) — 


— — — 20.1(1.8) 22.6(1.8) 17.6(2.1) 


1.7E4 18.6(10.8) 18.5(1.2) — 


— — — 18.5(1.2) 20.3(1.1) 16.8(1.4) 


3.5E4 19.7(4.6) 18.44(0.87) — 


— 19.15(0.70) 18.42(0.74) 18.39(0.90) 19.56(1.05) 17.34(0.70) 


7E4 18.42(0.43) 18.38(0.69) — 


— 18.82(0.61) 18.29(0.40) 18.33(0.69) 19.18(0.87) 17.64(0.69) 


1.3E5 18.41(0.26) 18.34(0.71) — 


— 18.72(0.55) 18.20(0.46) 18.28(0.72) 18.76(0.83) 17.78(0.80) 


2.7E5 18.27(0.21) 18.35(0.45) 18.47(1.03) 18.23(0.59) 18.55(0.42) 18.16(0.29) 18.29(0.45) 18.62(0.54) 18.09(0.46) 


5.5E5 18.26(0.13) 18.28(0.28) 18.25(0.49) 18.43(0.43) 18.44(0.32) 18.13(0.19) 18.20(0.29) 18.28(0.39) 18.25(0.26) 


1E6 18.23(0.13) 18.28(0.30) 18.23(0.30) 18.30(0.42) 18.32(0.26) 18.18(0.16) 18.21(0.31) 18.20(0.33) 18.25(0.31) 



TABLE III: Free energy difference estimates obtained for changing the charge of a Lennard- Jones particle from -e/2 to -|-e/2 
in a box of explicit water. Results are the averages from 16 independent simulations for various methods described in the text 
as a function of the number of dynamics steps used in the simulation. The standard deviation is shown in parentheses. For 
single-ensemble path sampling (SEPS and BSEPS) and Jarzynski methods (Jarz and BJarz), only the most efficient results 
are shown. The table shows that in the limit of long simulation times (10^ dynamics steps) all methods produce average AF 
estimates that roughly agree. The table also shows that AIM and BJarz approaches provide the most precise long-simulation 
estimate. 



that AIM has the abihty to produce more precise slope 
data compared to TI. 

Table IIIII shows AF estimates for the different ap- 
proaches. For all non-equilibrium approaches, only the 
most efficient data are shown. For SEPS and BSEPS 
the paths were composed of 500 A-steps (restricted to 
500 due to computer memory) with 80 dynamics steps 
at each value of A. For Jarz the paths were composed 
of 40 000 A-steps with one dynamics step at each value 
of A, and for BJarz, 20 000 A-steps with one dynamics 
step at each value of A were used. For all of these non- 



equilibrium data, the standard deviation of the work val- 
ues were aw ~ 0.8 kcal/mol ~ 1.3 fcsT, in agreement 
with previous studie o^^i^" , and with the growing data in 
this study. At least four different path lengths were at- 
tempted for each non-equilibrium method to determine 
the most efficient. 

Table Hill demonstrates that, for long simulation times, 
all methods produce roughly the same average AF es- 
timate. Also, the table shows that, given 10^ dynamics 
steps, AIM and BJarz methodologies provide the most 
precise free energy estimates. 
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FIG. 3: The slope of the free energy dF/dA as a function 
of A for a changing the charge of a Lennard- Jones particle 
in a box of explicit water from -e/2 to +e/2. Results for 
both TI and AIM methods are shown for 10*' dynamics steps. 
The data show the averages (data points) and standard devi- 
ations (error bars) from 16 independent simulations for each 
method. The errorbars are too small to resolve on the plot 
shown, however, it should be noted that the average uncer- 
tainty in the the slope for AIM is 0.38 kcal/mol and for TI 
is 1.05 kcal/mol, suggesting that AIM has the ability to pro- 
duce a more precise profile compared to TI. Thus, AIM is 
preferred over TI for computing the potential of mean force 
for this system. The smoothness of the profile also suggests 
that the switching function (7a of Eq. 1241 used in this report 
is adequate. 



Tables IIIII and IIVI show the difference between using 
equihbrium data in the "forward" (FEPF) and "reverse" 
(FEPR) directions. While, the results are similar for 10^ 
dynamics steps, it is clear that FEPF produces the de- 
sired results more rapidly than FEPR indicating that the 
configurational overlap is greater in the forward direc- 
tion. However, the FEPF data also tends to "overshoot" 
the correct value by a small margin which makes conver- 
gence of the FEPF estimate difficult to judge. 

For fast estimation of free energy differences, Table 
IIVI shows the number of dynamics steps needed by each 
method to obtain a free energy estimate within a spe- 
cific tolerance of AFiong sim (average of all estimates at 
10^ dynamics steps). Note that the number of dynamics 
steps needed for the SEPS and BSEPS methods are large 
due to the fact that many paths must be discarded for 
equilibration of the path ensemble. For all methods ex- 
cept AIM, the entries in Table Hvl were estimated using 
linear interpolation of the data in Table lllll From the 
data in the table, if the desired precision is less than 1.0 
kcal/mol, then all methods other than SEPS and BSEPS 
produce comparable results. However, if the desired pre- 
cision is less than 0.5 kcal/mol, then AIM and BJarz 
approaches are best. 

We conclude that, when charging a Lennard- Jones par- 
ticle in explicit solvent, the preferred methodology de- 
pends upon the type of estimate one wishes to gener- 
ate. If a very high quality estimate is desired, then AIM 



Method 


Within 1.0 kcal/mol 


Within 0.5 kcal/mol 


AIM 


52 000 


64 000 


TI 


27 500 


243 000 


SEPS 


291 000 


515 000 


BSEPS 


399 000 


487 000 


Jarz 


40 000 


180 000 


BJarz 


40 000 


69 000 


Benn 


29 000 


245 000 


FEPF 


43 000 


335 000 


FEPR 


26 000 


252 000 



TABLE IV: Number of dynamics steps necessary to be within 
a specified tolerance of the correct result A_Fiong sim = 18.24 
kcal/mol, average AF estimate at 10^ dynamics steps for all 
methods, for charging a Lennard- Jones particle in explicit 
solvent. The first column is the method used to obtain the 
estimate. The second column is the number of dynamics steps 
needed to estimate AF within 1.0 kcal/mol of AFiong sim with 
an uncertainty less than 1.0 kcal/mol. The third column is the 
number of dynamics steps needed to obtain an estimate within 
0.5 kcal/mol with an uncertainty less than 0.5 kcal/mol. 



is the best choice, closely followed by BJarz. If a very 
rapid estimate of AF, with an uncertainty of less than 
1.0 kcal/mol, is desired, then then comparable results 
are seen using all methodologies except for SEPS and 
BSEPS. If the A_F estimate is to be within 0.5 kcal/mol, 
then AIM and BJarz are the best choices. 

Finally, if the desired result is the potential of mean 
force, then AIM will generate a much smoother curve 
than TI. 



C. A second look at a two-dimensional model 

Because SEPS proved orders of magnitude more effi- 
cient than TI and Jarz in the study of a two-dimensional 
models, we return to that model in an effort to under- 
stand the decreased effectiveness of SEPS in the present 
study. Specifically, we use the model from Ref. 0, but 
now for a wide range of conformational sampling barrier 
heights (fixed A), and then compare SEPS to TI, as in 
our original study. Note, that we use the term "con- 
formational sampling barrier" to distinguish it from the 
barrier in A-space. 

Some alterations to our approach in Ref. |y were nec- 
essary to provide a fair comparison in the context of the 
present report. The results in Ref. 9 were obtained for 
very short paths, large perturbations of the shoot point, 
and a conformational sampling barrier height of 14.0 
ksT. For consistency with the present studies, SEPS 
results were generated with no perturbation of the shoot 
point, much longer paths, and for a range of conforma- 
tional sampling barrier heights. Both TI and SEPS simu- 
lations utilized Brownian dynamics to propagate the sys- 
tem. For SEPS, paths were generated as described in the 
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Barrier (fcsT) SEPS long SEPS short TI 



1.0 


60 000 


200 000 


15 300 


2.0 


120 000 


500 000 


35 700 


4.0 


400 000 


1 000 000 


204 000 


6.0 


1 400 000 


1 400 000 


1 020 000 


8.0 


8 000 000 


1 600 000 


5 100 000 


10.0 


40 000 000 


2 400 000 


20 400 000 


12.0 


80 000 000 


4 000 000 


76 500 000 


14.0 


200 000 000 


10 000 000 


204 000 000 



TABLE V: Number of dynamics steps necessary to be within 
0.5 ksT of the analytical result for AF with a 0.5 ksT or less 
standard deviation for the two-dimensional model in-. The 
first column is the barrier height of the potential energy sur- 
face in ksT units. The second and third columns are the 
number of dynamics steps using SEPS with, respectively, 200 
work values and 20 000 work values. The fourth column is 
the number of dynamics steps using TI with using 51 equally 
spaced values of A. For both TI and SEPS, half of the gener- 
ated data were thrown away for equilibration. 



present report (but with no velocity) , and accepted with 
the probability given in Eq. (|19|) . 

Results for the two-dimensional model using SEPS and 
TI are shown in Table IVl The free energy change is for 
switching between a single-well potential and a double- 
well potential with a conformational barrier height in 
ksT units given in the first column. The next three 
columns give the number of dynamics steps needed for 
the AF estimate to be within 0.5 fc^T of the correct value 
with 0.5 ksT or smaller standard deviation (estimated 
over at least 100 trials): the second and third columns 
are for SEPS where either 200 (long trajectories) or 20 
000 (short trajectories) work values were generated with 
50% of the work values discarded for equilibration, and 
the fourth column is TI using 51 evenly spaced values of 
A with 50% of the data at each value of A discarded for 
equilibration. 

Table^clearly shows that, for very low conformational 
barrier height, TI is much more efficient than SEPS, and 
that the most efficient SEPS is obtained using longer 
paths and thus fewer work values. For increasing con- 
formational barrier heights, SEPS using long paths and 
TI become comparable, while SEPS using short paths be- 
comes the most efficient. For the largest conformational 
barrier height tested in this study (14.0 fc^T), SEPS us- 
ing short paths is at least 20 times more efficient than 
either TI or SEPS using long paths. 

Since the results for growing and charging an ion in 
solvent showed that TI was more efficient than SEPS, we 
suggest that the free energy landscapes for the molecular 
systems used in this study have rather modest conforma- 
tional sampling barriers^Si^. 



VI. CONCLUSIONS 

We have carefully studied several computational free 
energy difference (AF) methods, comparing efficiency 
and precision. The test cases used for the compari- 
son were relative solvation energy calculations involv- 
ing either a large change in the Lennard- Jones size or 
in the charge of a particle in explicit solvent. Specifi- 
cally, we compared: adaptive integration (AIM)^; ther- 
modynamic integration (TI)ii; path sampling of non- 
equilibrium work values using both a Jarzynski uni- 
directional formalism (SEPS)^, and a Bennett-like bi- 
directional formalism (BSEPS); Jarzynski (Jarz)^^ and 
Bennett (BJarz)^^'^^ averaging of non-equilibrium work 
values; equilibrium Bennett (Benn)i^; and free energy 
perturbation (forward, FEPF and reverse FEPR)— . 

AIM^ was found to provide very high quality, precise 
estimates, given long simulation times (10^ total dynam- 
ics steps in this study), and also allowed very rapid es- 
timation of AF. In addition, AIM provided smooth free 
energy profiles (and thus smooth potential of mean force 
curves) as compared to TI; see Figs. ^ and 13 Clearly, 
AIM was the best all-around choice for the systems stud- 
ied here. 

B JaraW was also found to perform very well, with long- 
simulation results that were second only to AIM. How- 
ever, it should be noted that the data shown in this study 
are for the most efficient path lengths only. To deter- 
mine the optimal path length, many simulations were 
performed, adding to the overall cost of the method. 
Also, our results showed that using bi-directional data 
(BJarz) produced considerably more precise results than 
using uni-directional data (Jarz). 

The SEPS method is shown to provide accurate free 
energy estimates from "fast-growth" uni-directional non- 
equilibrium work values. Specifically, in cases where the 
standard deviation of the work values is much greater 
than ksT {aw ^ ksT), the SEPS method can effec- 
tively shift the work values to allow for more accurate AF 
estimation than is possible using ordinary Jarzynski av- 
eraging. Interestingly, using bi-directional data (BSEPS) 
did not increase the precision of the AF estimate, and 
perhaps made it somewhat worse. 

We also find, in agreement with previous studiesi^iSS, 
that the greatest efficiency for the Jarz approach is when 
aw ~ 1 ksT. For the first time, we also show that SEPS 
is also most efficient when aw ~ 1 ksT, for the systems 
studies in this report. 

We have also suggested an explanation — with poten- 
tially quite interesting consequences — for the decreased 
effectiveness of SEPS in molecular systems. By re- 
examining the two-dimensional model used in our first 
SEPS paper^, we find that SEPS can indeed be much 
more more efficient than TI, but only when the confor- 
mational sampling barrier is very high (S> ksT). This 
suggests that the configurational sampling barriers en- 
countered in the molecular systems studied here are fairly 
modest, counter to our own expectations. A key ques- 
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tion is thus raised: How high are conformational sam- 
pUng barriers encountered in free energy calculations of 
"practical interest?" See also Refs. li^lHfl 

We remind the reader that the results of this study 
are valid only for the types of AF calculations we 
considered — namely, growing and charging a Lennard- 
Jones particle in explicit solvent. When large confor- 
mational changes are important, such as for binding 
affinities, the results could be significantly different — 
particularly if large conformational sampling barriers are 
present. 
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